  
/*********************************************************************
*                                                                    *
* F -- evaluates the function f(y,t), lots of physical imput on it   *
*                                                                    *
* Parameters:                                                        *
*       fields -- pointer to field_vector from where to extract (y,t)*
*       derivs -- pointer to field_vector where to put derivatives   *
*     wave_par -- pointer to parameter struct                        *
*                                                                    *
* Returns: nothing                                                   *
*                                                                    *
*********************************************************************/
#include "first_macro_2d.h"  /* Where global parameters are defined */
#include "structs_2d.h"      /* Where structures are defined */
#include "derivs_2d.h"       /* Where derivatives functions are defined */
 
#include "gen_2d.h"

#define GRD_IND [IND_2(grid_ind1,grid_ind2)]


void heat_2d_eq(struct GRID_PAR *grid_ptr,
		struct field_array *fields_ptr, 
		struct field_array *derivs_ptr,
		struct FUNCTION_PAR *function_par) {


    int ni_1 = (*grid_ptr).start_grid_1; 
    int nf_1 = (*grid_ptr).final_grid_1;  
    int ni_2 = (*grid_ptr).start_grid_2;  
    int nf_2 = (*grid_ptr).final_grid_2; 
    int n_gridpts_1 = nf_1 - ni_1;
    int n_gridpts_2 = nf_2 - ni_2;
    int grid = (*grid_ptr).grid;
    int coorx = 1;    
    int coory = 2;


    FLOAT xi = (*grid_ptr).initial_x;
    FLOAT xf = (*grid_ptr).final_x;
    FLOAT yi = (*grid_ptr).initial_y;
    FLOAT yf = (*grid_ptr).final_y;

//    FLOAT x, y, m_t_nn, nn_x, nn_y, nn_z;

    FLOAT one_dN1 = (xf-xi)/(FLOAT)(nf_1-ni_1-1);
    FLOAT one_dN2 = (yf-yi)/(FLOAT)(nf_2-ni_2-1);

    FLOAT sigma = (*function_par).sigma; /* For dissipation */

    FLOAT mm = (*function_par).mm;

    FLOAT delta = (*function_par).delta; // used in penalty

    /* normals */

    FLOAT nx_x0 = (*grid_ptr).nx[0][0];
    FLOAT ny_x0 = (*grid_ptr).ny[0][0];
    
    FLOAT nx_x1 = (*grid_ptr).nx[0][1];
    FLOAT ny_x1 = (*grid_ptr).ny[0][1];
    
    FLOAT nx_y0 = (*grid_ptr).nx[1][0];
    FLOAT ny_y0 = (*grid_ptr).ny[1][0];
  
    FLOAT nx_y1 = (*grid_ptr).nx[1][1];
    FLOAT ny_y1 = (*grid_ptr).ny[1][1];






    /* Corner values */





    FLOAT nx;
    FLOAT ny;

    FLOAT m_x=0.0;
    FLOAT dx_m_x=0.0;
    FLOAT dy_m_x=0.0;
    FLOAT psi_x=0.;
    FLOAT m_y=0.0;
    FLOAT dx_m_y=0.0;
    FLOAT dy_m_y=0.0;
    FLOAT psi_y=0.;
    FLOAT m_z=0.0;
    FLOAT dx_m_z=0.0;
    FLOAT dy_m_z=0.0;
    FLOAT psi_z=0.;
    FLOAT d2_x, d2_y, d2_z, d2, md_xm, md_ym, mmd2,norm_m, m_t_dot_m ;

    FLOAT m_x_n=0.0;
    FLOAT psi_x_n=0.0;
    FLOAT m_y_n=0.0;
    FLOAT psi_y_n=0.0;
    FLOAT m_z_n=0.0;
    FLOAT psi_z_n=0.0;

    FLOAT m_x_t=0.0;
    FLOAT m_y_t=0.0;
    FLOAT m_z_t=0.0;

    FLOAT norm_n, detg_1_2;

    FLOAT g_up_xx, g_up_yy, g_up_xy, g_up_yx;


    FLOAT dm_x, dm_y, dm_z, dpsi_x, dpsi_y, dpsi_z, Tau;

    FLOAT  factor, o_dx;

  char macro_value_strg[100];



    GET_MACRO_VALUE(DERIV);
    if (strcmp(macro_value_strg,"derivS_2d")==0) {
      factor = 2.0;
    }
    else if (strcmp(macro_value_strg,"derivQ_2d")==0) {
      factor = 48.0/17.0;
    }
    else if (strcmp(macro_value_strg,"derivQ_3_2d")==0) {
      factor = 11.0/3.0;
    }
    else if (strcmp(macro_value_strg,"deriv_strand_third_order_boundaries_sixth_interior_2d")==0) {
      factor = 43200.0/13649.0;
    }
    else if (strcmp(macro_value_strg,"deriv_strand_fourth_order_boundaries_eight_interior_2d")==0) {
      factor = 5080320.0/1498139.0;
    }
    else {
      factor=2.0; printf("check factor por penalty!!!%s!!!",macro_value_strg);
    }



    
    o_dx = factor*(FLOAT)(nf_1 - ni_1 -1)/(xf-xi);
    Tau= -(1.0/8.0)*o_dx*((1.0+delta)*(1.0+delta)+delta*delta);
    FLOAT A = 1.0;
















  /* first the time */

  (*derivs_ptr).time = (*fields_ptr).time;



// PROJECTION INTO THE UNIT SPHERE

#ifdef PROJECTION

{register int grid_ind1, grid_ind2;
 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
   for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){


       	m_x = (*fields_ptr).u[M_A]GRD_IND;
       	m_y = (*fields_ptr).u[M_B]GRD_IND;
       	m_z = (*fields_ptr).u[M_C]GRD_IND;

	norm_m = sqrt(m_x*m_x + m_y*m_y + m_z*m_z);
        
	(*fields_ptr).u[M_A]GRD_IND = (*fields_ptr).u[M_A]GRD_IND/norm_m;
	(*fields_ptr).u[M_B]GRD_IND = (*fields_ptr).u[M_B]GRD_IND/norm_m;
	(*fields_ptr).u[M_C]GRD_IND = (*fields_ptr).u[M_C]GRD_IND/norm_m;
	

	 }
     }
  }
 
#endif


  /* MAKE swap fields AND take derivatives */

#ifdef DEBUG
 printf("in eqn. before taking derivatives \n");
#endif





DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_A],
         (struct field *)&dfields.du[D_X_M_A],&coorx);
DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_A],
         (struct field *)&dfields.du[D_Y_M_A],&coory);

DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_B],
         (struct field *)&dfields.du[D_X_M_B],&coorx);
DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_B],
         (struct field *)&dfields.du[D_Y_M_B],&coory);

DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_C],
         (struct field *)&dfields.du[D_X_M_C],&coorx);
DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_C],
         (struct field *)&dfields.du[D_Y_M_C],&coory);


//     M_A    

 {register int grid_ind1, grid_ind2;
 
 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
   for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){

         g_up_xx = auxfields.u_aux[G_UP_XX]GRD_IND;
         g_up_yy = auxfields.u_aux[G_UP_YY]GRD_IND;
         g_up_xy = auxfields.u_aux[G_UP_XY]GRD_IND;

	 detg_1_2 = 1./sqrt(g_up_xx*g_up_yy - g_up_xy*g_up_xy);

	 auxfields.u_aux[SWAP_X_X]GRD_IND = 
	 (g_up_xx * dfields.du[D_X_M_A]GRD_IND + g_up_xy * dfields.du[D_Y_M_A]GRD_IND)*detg_1_2;
	 auxfields.u_aux[SWAP_X_Y]GRD_IND = 
	 (g_up_yy * dfields.du[D_Y_M_A]GRD_IND + g_up_xy * dfields.du[D_X_M_A]GRD_IND)*detg_1_2; 


	 auxfields.u_aux[SWAP_Y_X]GRD_IND = 
	 (g_up_xx * dfields.du[D_X_M_B]GRD_IND + g_up_xy * dfields.du[D_Y_M_B]GRD_IND)*detg_1_2;
	 auxfields.u_aux[SWAP_Y_Y]GRD_IND = 
	 (g_up_yy * dfields.du[D_Y_M_B]GRD_IND + g_up_xy * dfields.du[D_X_M_B]GRD_IND)*detg_1_2; 


	 auxfields.u_aux[SWAP_Z_X]GRD_IND = 
	 (g_up_xx * dfields.du[D_X_M_C]GRD_IND + g_up_xy * dfields.du[D_Y_M_C]GRD_IND)*detg_1_2;
	 auxfields.u_aux[SWAP_Z_Y]GRD_IND = 
	 (g_up_yy * dfields.du[D_Y_M_C]GRD_IND + g_up_xy * dfields.du[D_X_M_C]GRD_IND)*detg_1_2;


     }
   }
 }


{register int  grid_ind1, grid_ind2; 

     for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){

	 grid_ind1 = ni_1;
 	 sfield_send.s_x_0_u[M_A][grid_ind2] = (*fields_ptr).u[M_A]GRD_IND;
 	 sfield_send.s_x_0_u[M_B][grid_ind2] = (*fields_ptr).u[M_B]GRD_IND;
 	 sfield_send.s_x_0_u[M_C][grid_ind2] = (*fields_ptr).u[M_C]GRD_IND;

 	 sfield_send.s_x_0_u[PSI_X][grid_ind2] = -auxfields.u_aux[SWAP_X_X]GRD_IND;
 	 sfield_send.s_x_0_u[PSI_Y][grid_ind2] = -auxfields.u_aux[SWAP_Y_X]GRD_IND;
 	 sfield_send.s_x_0_u[PSI_Z][grid_ind2] = -auxfields.u_aux[SWAP_Z_X]GRD_IND;

	 grid_ind1= nf_1-1;
	 sfield_send.s_x_1_u[M_A][grid_ind2] = (*fields_ptr).u[M_A]GRD_IND;	 sfield_send.s_x_1_u[M_B][grid_ind2] = (*fields_ptr).u[M_B]GRD_IND;
	 sfield_send.s_x_1_u[M_C][grid_ind2] = (*fields_ptr).u[M_C]GRD_IND;

 	 sfield_send.s_x_1_u[PSI_X][grid_ind2] = auxfields.u_aux[SWAP_X_X]GRD_IND; 	 sfield_send.s_x_1_u[PSI_Y][grid_ind2] = auxfields.u_aux[SWAP_Y_X]GRD_IND;
 	 sfield_send.s_x_1_u[PSI_Z][grid_ind2] = auxfields.u_aux[SWAP_Z_X]GRD_IND;

     }
 
      
     for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
	 grid_ind2 = ni_2;
	 sfield_send.s_y_0_u[M_A][grid_ind1] = (*fields_ptr).u[M_A]GRD_IND;	 sfield_send.s_y_0_u[M_B][grid_ind1] = (*fields_ptr).u[M_B]GRD_IND;
	 sfield_send.s_y_0_u[M_C][grid_ind1] = (*fields_ptr).u[M_C]GRD_IND;

	 sfield_send.s_y_0_u[PSI_X][grid_ind1] = -auxfields.u_aux[SWAP_X_Y]GRD_IND;	 sfield_send.s_y_0_u[PSI_Y][grid_ind1] = -auxfields.u_aux[SWAP_Y_Y]GRD_IND;
	 sfield_send.s_y_0_u[PSI_Z][grid_ind1] = -auxfields.u_aux[SWAP_Z_Y]GRD_IND;


	 grid_ind2 = nf_2-1;
	 sfield_send.s_y_1_u[M_A][grid_ind1] = (*fields_ptr).u[M_A]GRD_IND;	 sfield_send.s_y_1_u[M_B][grid_ind1] = (*fields_ptr).u[M_B]GRD_IND;
	 sfield_send.s_y_1_u[M_C][grid_ind1] = (*fields_ptr).u[M_C]GRD_IND;

	 sfield_send.s_y_1_u[PSI_X][grid_ind1] = auxfields.u_aux[SWAP_X_Y]GRD_IND;
	 sfield_send.s_y_1_u[PSI_Y][grid_ind1] = auxfields.u_aux[SWAP_Y_Y]GRD_IND;
	 sfield_send.s_y_1_u[PSI_Z][grid_ind1] = auxfields.u_aux[SWAP_Z_Y]GRD_IND;

     }
 }






DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_X_X],
         (struct field *)&dfields.du[DD_X_M_A],&coorx);

DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_X_Y],
         (struct field *)&dfields.du[DD_Y_M_A],&coory);

DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_Y_X],
         (struct field *)&dfields.du[DD_X_M_B],&coorx);

DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_Y_Y],
         (struct field *)&dfields.du[DD_Y_M_B],&coory);

DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_Z_X],
         (struct field *)&dfields.du[DD_X_M_C],&coorx);

DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_Z_Y],
         (struct field *)&dfields.du[DD_Y_M_C],&coory);






#ifdef DISSIPATION

DISS(grid_ptr,(struct field *)(*fields_ptr).u[M_A],
         (struct field *)&dfields.du[DISS_M_A]);

DISS(grid_ptr,(struct field *)(*fields_ptr).u[M_B],
         (struct field *)&dfields.du[DISS_M_B]);

DISS(grid_ptr,(struct field *)(*fields_ptr).u[M_C],
         (struct field *)&dfields.du[DISS_M_C]);
#endif

#ifdef DEBUG
   printf("in eqn. after taking derivatives \n");
#endif





    /* ---------------------------------- PMI part --------------------------------------------- */

#include"MPI_SR.c"
 /* --------------------------------------------------------------------------------------------- */


/*************************************************************************************************/
/*************************************************************************************************/

/* inner points */

/*************************************************************************************************/

 {register int grid_ind1, grid_ind2;
 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
   for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){


#include"heat_2d_eq_eq.c"


	     (*derivs_ptr).u[M_A]GRD_IND =  m_x_t;
	     (*derivs_ptr).u[M_B]GRD_IND =  m_y_t;
	     (*derivs_ptr).u[M_C]GRD_IND =  m_z_t;


	 }
     }
  }
 

/* --------------------------------------------- MPI -------------------------------------- */
/* ---------------------------------- check fields have been received and sent ------------ */

#include"MPI_W.c"

/* ------------------------------------------------------------------------------------------------------*/
/* ------------------------------------SWAP order of some faces------------------------------------------*/
/* ------------------------------------------------------------------------------------------------------*/

#include"MPI_S2_inv.c"

/* ------------------------------------------------------------------------------------------------------*/
/* --------------------------------SWAP direction of some fields-----------------------------------------*/
/* ------------------------------------------------------------------------------------------------------*/
 
//#include"wave_2d_eq_jaco.c"

/*-------------------------------------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------------------------------------*/

 /*----------------------------- BOUNDARY CONDITIONS --------------------------*/

  /* Boundary conditions at x=xi and x=xf */

	 /************************* x = xi ******************************/

{register int grid_ind1;
 register int grid_ind2;
 
 grid_ind1 = ni_1;

       nx=nx_x0;
       ny=ny_x0;

  
 for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){

 m_x_n = sfield_recv.s_x_0_u[M_A][grid_ind2];
 psi_x_n = sfield_recv.s_x_0_u[PSI_X][grid_ind2];
 m_y_n = sfield_recv.s_x_0_u[M_B][grid_ind2];
 psi_y_n = sfield_recv.s_x_0_u[PSI_Y][grid_ind2];
 m_z_n = sfield_recv.s_x_0_u[M_C][grid_ind2];
 psi_z_n = sfield_recv.s_x_0_u[PSI_Z][grid_ind2];

#include"heat_2d_eq_side.c"

     }
 }

  
	 /************************* x = xf ******************************/

{register int grid_ind1;
 register int grid_ind2;


      nx=nx_x1;
      ny=ny_x1;

      grid_ind1 = nf_1-1;

 for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){

 m_x_n = sfield_recv.s_x_1_u[M_A][grid_ind2];
 psi_x_n = sfield_recv.s_x_1_u[PSI_X][grid_ind2];
 m_y_n = sfield_recv.s_x_1_u[M_B][grid_ind2];
 psi_y_n = sfield_recv.s_x_1_u[PSI_Y][grid_ind2];
 m_z_n = sfield_recv.s_x_1_u[M_C][grid_ind2];
 psi_z_n = sfield_recv.s_x_1_u[PSI_Z][grid_ind2];

#include"heat_2d_eq_side.c"

  }
  }
 


  /* Boundary conditions at y = yi and y = yf */

	 /************************* y = yi ******************************/

 {register int grid_ind1;
 register int grid_ind2;

       nx=nx_y0; 
       ny=ny_y0; 

       grid_ind2 = ni_2;

 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){

 m_x_n = sfield_recv.s_y_0_u[M_A][grid_ind1];
 psi_x_n = sfield_recv.s_y_0_u[PSI_X][grid_ind1];
 m_y_n = sfield_recv.s_y_0_u[M_B][grid_ind1];
 psi_y_n = sfield_recv.s_y_0_u[PSI_Y][grid_ind1];
 m_z_n = sfield_recv.s_y_0_u[M_C][grid_ind1];
 psi_z_n = sfield_recv.s_y_0_u[PSI_Z][grid_ind1];

#include"heat_2d_eq_side.c"

 }
 }


	 /************************* y = yf ******************************/

 {register int grid_ind1;
 register int grid_ind2;

      nx=nx_y1;
      ny=ny_y1;

	 grid_ind2 = nf_2-1;

 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){

 m_x_n = sfield_recv.s_y_1_u[M_A][grid_ind1];
 psi_x_n = sfield_recv.s_y_1_u[PSI_X][grid_ind1];
 m_y_n = sfield_recv.s_y_1_u[M_B][grid_ind1];
 psi_y_n = sfield_recv.s_y_1_u[PSI_Y][grid_ind1];
 m_z_n = sfield_recv.s_y_1_u[M_C][grid_ind1];
 psi_z_n = sfield_recv.s_y_1_u[PSI_Z][grid_ind1];

#include"heat_2d_eq_side.c"
 
 }
  }
 
  

   
  
/***********************************************************************/

 /*-------------------------------- Corner Conditions -------------------------------*/

/***********************************************************************/

 /* WARNING corner conditions have been elliminated by multiplying by zero. */

  /* Corner x=0,y=0 */
//  {register int grid_ind1, grid_ind2, grid_ind3;
//  FLOAT SP_int_x, SM_int_x, ONEoverSPMSM_int_x, rN_int_x, V_p_int_x;
//  FLOAT SP_int_y, SM_int_y, ONEoverSPMSM_int_y, rN_int_y, V_p_int_y;
// 
//  grid_ind1 = ni_1;
//  grid_ind2 = ni_2;
// 
//  nx=nx_x0;
//  ny=ny_y0;
// 
//  phi_n = sfield_recv.s_x_0_u[PHI][grid_ind2];
//  psi_n = sfield_recv.s_x_0_u[PSI][grid_ind2];
// 
//  phi_t_n_y = sfield_recv.s_y_0_u[PHI][grid_ind1];
//  phi_x_n_x = sfield_recv.s_y_0_u[PSI][grid_ind1];
//  phi_y_n_y = sfield_recv.s_y_0_u[PHI_Y][grid_ind1];
// 
// #include"wave_2d_eq_corners.c"
// 
//  }



  /* Corner x=1,y=0 */
//  {register int grid_ind1 = nf_1-1;
//  register int grid_ind2 = ni_2;
//  nx= nx_x1; /* normal value nx=1 */
//  ny= ny_y0;   /* normal value ny=0 */
// 
//  phi_t_n_x = sfield_recv.s_x_1_u[PHI_T][grid_ind2];
//  phi_x_n_x = sfield_recv.s_x_1_u[PHI_X][grid_ind2];
//  phi_y_n_x = sfield_recv.s_x_1_u[PHI_Y][grid_ind2];
// 
//  phi_t_n_y = sfield_recv.s_y_0_u[PHI_T][grid_ind1];
//  phi_x_n_y = sfield_recv.s_y_0_u[PHI_X][grid_ind1];
//  phi_y_n_y = sfield_recv.s_y_0_u[PHI_Y][grid_ind1];
// 
// 
// #include"wave_2d_eq_corners.c"
// 
// }



  /* Corner x=0,y=1,z=0 */
//  {register int grid_ind1 = ni_1;
//  register int grid_ind2 = nf_2-1;
//  nx= nx_x0; /* normal value nx=0 */
//   ny= ny_y1;   /* normal value ny=0 */
// 
//  phi_t_n_x = sfield_recv.s_x_0_u[PHI_T][grid_ind2];
//  phi_x_n_x = sfield_recv.s_x_0_u[PHI_X][grid_ind2];
//  phi_y_n_x = sfield_recv.s_x_0_u[PHI_Y][grid_ind2];
// 
//  phi_t_n_y = sfield_recv.s_y_1_u[PHI_T][grid_ind1];
//  phi_x_n_y = sfield_recv.s_y_1_u[PHI_X][grid_ind1];
//  phi_y_n_y = sfield_recv.s_y_1_u[PHI_Y][grid_ind1];
// 
// #include"wave_2d_eq_corners.c"
// 
//  }


  /* Corner x=1,y=1,z=0 */
// {register int grid_ind1 = nf_1-1;
//  register int grid_ind2 = nf_2-1;
//   nx= nx_x1; /* normal value nx=-1 */
//   ny= ny_y1;   /* normal value ny=0 */
// 
//  phi_t_n_x = sfield_recv.s_x_1_u[PHI_T][grid_ind2];
//  phi_x_n_x = sfield_recv.s_x_1_u[PHI_X][grid_ind2];
//  phi_y_n_x = sfield_recv.s_x_1_u[PHI_Y][grid_ind2];
// 
//  phi_t_n_y = sfield_recv.s_y_1_u[PHI_T][grid_ind1];
//  phi_x_n_y = sfield_recv.s_y_1_u[PHI_X][grid_ind1];
//  phi_y_n_y = sfield_recv.s_y_1_u[PHI_Y][grid_ind1];
// 
// #include"wave_2d_eq_corners.c"
// 
//  }

// fflush(stdout);
//  exit(0);

/*  printf("At the end of eq. in n%d \n", grid); */
}
/*************************************************************************************/

 




